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1.  Definitions 

Let  A be  a rectangular  real  (mxn)  matrix,  and  X be  a 
(nxm)  real  matrix  satisfying  the  following  conditions  (t  de- 
notes transpose): 


XAX  = X 


(XAf  = XA  (4) 

If  X satisfies  only  (1)  it  is  called  a generalized  in- 

91 

verse,  and  will  be  denoted  by  A = X.  [1];  when  X satisfies 

g2 

/ 1 \ /o\  i v _ n * - i • j • 


(1)  and  (2)  we  call  X 


a reflexive  generalized  inverse; 


when  X satisfies  (1),  (2),  (3)  and  (4)  it  is  denoted  by  A+  and 
is  called  the  Moore-Penrose  pseudoinverse. 

The  partitioning  of  a matrix  A into  submatrices  Aj^  and 
A2  along  the  column  (or  row)  direction  will  be  denoted  by 


- AilA; 


( M 

V A2/ 


The  computation  of  the  pseudoinverse  of  a large  matrix 
is  very  important  for  several  applications  in  image  processing 
and  filtering  [2] . Since  the  matrices  involved  in  these  prob- 
lems are  very  large,  it  becomes  necessary  to  use  block  par- 
titioning schemes  for  this  purpose.  Unfortunately,  however, 
the  various  methods  used  for  block  partitioning  have  several 


theoretical  limitations,  if  one  wants  to  set  up  a recursive 


Then 


For  these  expressions  to  be  valid,  the  following  con 


ditions  are  necessary: 


*2 
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(i)  A should  be  positive  semi -definite  for  using  (6) 

(ii)  A should  be  positive  semi -definite  and  Rank  (A)  * 
Rank  (A^)  + Rank  (A22)  for  using  (7)  . 

Therefore,  if  we  want  to  compute  A+  recursively,  by  a 
successive  partitioning  procedure  involving  smaller  matrices, 
then  at  any  stage  the  intermediate  partitions  or  block  matrices 
should  also  satisfy  the  above  conditions.  It  is  clear  that 
while  it  is  possible  to  make  A positive  definite  by  starting 
with  AAfc  and  computing  A+  using  the  formula 

A+  « At(AAt*AAt)92AAt 

it  cannot  always  be  guaranteed  that  every  further  su  - 
partition  would  satisfy  the  above  conditions.  Therefore,  in 
order  to  recursively  partition  A through  several  stages,  we 
must  insure  that  at  least  the  semi -definiteness  condition  is 
satisfied  by  suitable  multiplication.  This  makes  the  partition 
ing  procedure  very  complex. 

In  order  to  obviate  this  difficulty,  a recursive  bi- 
section-transposition algorithm  is  suggested  below  based  on 
the  available  results  on  g-inverses  [1] . This  algorithm  does 
not  depend  upon  the  conditions  (i)  and  (ii)  since  at  any  stage 
either  a row  or  column  partitioning  of  the  matrix  is  carried 
out  (and  not  both) . 


When  B_  is  a 


3.  Bisection-Transposition  Algorithm 


a.  Principle 


Let  the  given  matrix  Ar-1  (mxn)  be  partitioned  as 


A , = [A  B ] 
r— 1 r1  r 


where  Ar(mxs)  and  Bf  (mx(n-s)Ks<n) . Then  it  is  easily  proved 


that  [1] 


Remarks 


(i)  When  Ar  is  non-singular,  we  get  Cr  = 0 and  Ar  = A 


column  vector,  Cr  is  also  a on  ■>  vector  and 


hence  CrCr  ■ 1 and  hence  Dr 


Therefore,  we  get 


where 


(11) 


(12) 


(13) 


ArBrCr 


Equation  (16)  corresponds  to  the  discrete  Kalman 
filtering  equations.  (Kalman  filtering  is  used  ex- 
tensively to  estimate  the  internal  state  variable 
of  a linear  system  based  on  noisy  measurements  of 
output  variables.) 

ii)  In  (12)  and  (13)  certain  matrices  need  not  be  con- 
formable; this  is  taken  care  of  by  appending  the 
required  zero  rows  and  columns  to  the  smaller  matrices 


b.  Algorithm 

Let  Aq  be  the  given  (mxm)  matrix  and  let  (pxp)  be 
the  order  to  which  we  finally  desire  to  reduce  this  matrix. 
Then  the  number  of  partitions  involved  is  r = 2 log,  (-)  ; the 

number  of  transpositions  required  is  (r-1) . The  algorithm 
uses  any  one  of  the  standard  procedures  for  computing 

A+  [4,  5]. 

In  this  algorithm  -*•  indicates  the  partitioning  operation 
and  the  use  of  (10)  to  obtain  the  pseudoinverse  of  the  pre- 
ceding larger  matrix. 


Step  Is  Set  i - 1;  Aq  [Ai|Bil 
Step  2s  A*  - IAi+ilBi+il 


" vs^.tvjuwi  n. 


yr  v. - vA.«r-r.>- v«*-v  • s - 


Step  3s  Set  i = i+1 

Step  4:  Is  i * r;  if  yes  go  to  Step  2;  otherwise  go  go  Step 


Step  5 
Step  6 
Step  7 
Step  8 


Set  j = r;  compute  At  . 

(Aj-l)+  ■"  <Aj  lBj>  • 
sot  j = j-1 

Is  j = 0?  If  yes,  go  to  Step  9;  otherwise  go  to  Step 


Step  9:  Result  = AQ.  Stop. 


Note: 


In  partitioning,  it  is  not  required  to  carry  out  exact 


bisection  of  the  matrix  A^;  however,  it  is  convenient  to  have 
exact  bisection,  and  to  have  the  order  of  AQ  be  a power  of  2. 


110  1 
1110 
0 111 
10  11 


Aq  is  not  positive  semi -definite  and  r(Ai;L)  + r(A22)  ? r^A) 


If  we  take  A,  = 11  ; A. 

11  12  10 


A2l“  10  ; A22 


El- 


Therefore,  the  classi- 


cal four  partitioning  scheme  will  fail. 

We  now  use  the  successive  bisection-transposition  al- 
gorithm; here  m=  4,  p = 2,  r = 2 . 


Take 


; B, 


fnoil  rn* 

Lllioj  f A2  = Lll. 


[::] 


IK] 


(I  - a2a2)  B: 


a 


,.  T ■ * ■ :ih 


. 


c+  = i r-i  -i] 

^2  2 L 1 -lj 


(I  - C2C+2)  B^(A+)t 


H2H2t  mk[l\] 


(I  - CJC2)Q-J  B*(A*)aJ(I 


T [l  l] 


B2C2  * “ 10  ( 1 1 1 


_ i ri  i -2  3] 

5 LI  1 3 -2 J 


Using  A^,  A^  and  we  get 


(I  - AlAl)  B, 


-4  1 

11-4 
5 3 3 

3 3 


+ _ 1 [-2  1 1 11 
C1  ‘ 1 L 1 -2  1 lj 


CjC^  =*  0;  D x 


, 11-21 
± 1 1 1-2 

J -2  1 1 1 

1-211 


. jt 


Concluding  Remarks 

(i)  The  above  algorithm  is  useful  when  one  deals  with 
matrices  having  certain  special  block  features. 

(ii)  A fast  algorithm  for  transposition  is  available  in 

[6]. 

iii)  The  numerical  stability  of  the  above  procedure  is 
not  established. 

It  is  well-known  that  the  numerical  algorithms 
for  computing  pseudoinverses  lack  numerical  stability 
and  accordingly,  error-free  exact  calculations  are 
desirable  [5];  see  also  [7]. 


(iv)  The  computational  complexity  of  the  above  algorithm 
is  to  be  established  [8] . 
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